function riMain(inDir,outDir)

%inDir='C:\Users\matthew\Desktop\2012\calc\results\';
%retr_dir='retrieval_theta';
%proc_dir='procedural_theta';
%stat_dir='C:\Users\matthew\Desktop\2012\calc\results\stats\theta\';
nSurr=100;
files=dir(fullfile(inDir,'*alpha*.mat'));
for i=1:length(files)
    filename=files(i).name;
    load(fullfile(inDir,filename),'pliXSubjs');
    %thSubjData=thresholdData(pliXSubjs,.75);
    [smallworld,modularity]=graphAnalysis(pliXSubjs,nSurr);
    outfile=strtok(filename,'.');
    outfile=strcat('sw_mod_',outfile);
    save(fullfile(outDir,outfile),'smallworld','modularity');
end


%global_stat(smallworld_one,smallworld_two,modularity_one,modularity_two);


function [smallworld,modularity]=graphAnalysis(subjData,nSurr)

nSubjs=size(subjData,1);
nNodes=size(subjData,2);
%nSurr=50;

smallworld.C=zeros(nSubjs,1);
smallworld.lambda=zeros(nSubjs,1);
smallworld.efficiency=zeros(nSubjs,1);

smallworld.Csurr=zeros(nSubjs,1);
smallworld.lambdasurr=zeros(nSubjs,1);
smallworld.efficiencysurr=zeros(nSubjs,1);

modularity.Q=zeros(nSubjs,1);
modularity.Ci=zeros(nSubjs,nNodes);
modularity.PC=zeros(nSubjs,nNodes);
modularity.Z=zeros(nSubjs,nNodes);

modularity.Qsurr=zeros(nSubjs,1);
modularity.PCsurr=zeros(nSubjs,nNodes);
modularity.Zsurr=zeros(nSubjs,nNodes);

for j=1:nSubjs

    m=squeeze(subjData(j,:,:));
   
    [smallworld.C(j),smallworld.lambda(j),smallworld.efficiency(j)]=smallworld_indexes(m);

       
    [Ci,Q]=modularity_und(m);
    modularity.Q(j)=Q;
    modularity.Ci(j,:)=Ci';
    pc=participation_coef(m,Ci);
    modularity.PC(j,:)=pc';
    z=module_degree_zscore(m,Ci);
    modularity.Z(j,:)=z';
    
    
   [ smallworld.Csurr(j), smallworld.lambdasurr(j),smallworld.efficiencysurr(j), ...
       modularity.Qsurr(j),modularity.PCsurr(j,:),modularity.Zsurr(j,:)]=surrogateNet(m,nSurr);
 
end

   
function [surrC,surrLambda,surrEfficiency,surrQ,surrPC,surrZ]=surrogateNet(m,nSurr)

nCh=size(m,1);

cV=zeros(nSurr,1);
lambdaV=zeros(nSurr,1);
efficiencyV=zeros(nSurr,1);

Qv=zeros(nSurr,1);
PC=zeros(nSurr,nCh);
Z=zeros(nSurr,nCh);

for i=1:nSurr
    
    rw=null_model_und_sign(m,1);
    %rw=randmio_und_connected(m,1);

    [cV(i),lambdaV(i),efficiencyV(i)]=smallworld_indexes(rw);
    [Ci,Qv(i)]=modularity_und(rw);
    pc=participation_coef(rw,Ci);
    PC(i,:)=pc';
    z=module_degree_zscore(rw,Ci);
    Z(i,:)=z';
end

surrC=mean(cV,1);
surrLambda=mean(lambdaV,1);
surrEfficiency=mean(efficiencyV,1);
surrQ=mean(Qv,1);
surrPC=mean(PC,1);
surrZ=mean(Z,1);
surrPC=surrPC';
surrZ=surrZ';

function [c,lambda,efficiency]=smallworld_indexes(conMatrix)

c = clustering_coef_wu(conMatrix); %average cluster coefficient for entire net   
c=nanmean(c);
conMatrix2d=1./conMatrix;             %from correlation values to lengths
D=distance_wei(conMatrix2d);       %distance matrix of min length path
[lambda,efficiency] = charpath(D);
%[lambda,efficiency]=Mycharpath(D); %charateristic path (lambda) e efficiency

function [areas,arealabels]=areasAndLabels()

[arealabels]=getAreaLabels();
[areas]=getAreas();
function [arealabels]=getAreaLabels()
arealabels={'laf';'lf';'lfc';'lc';'lcp';'lp';'lpo';'lt';'raf';'rf';'rfc';'rc';'rcp';'rp';'rpo';'rt'};
function [areas]=getAreas()

laf=[22 25 26];
lf=[30 31 32 33];
lfc=[1 2 3];
lc=[8 9 10];
lcp=[15 16 17];
lp=[47 48 49 50];
lpo=[56 57 61];
lt=[39 41 43 45];

raf=[24 28 29];
rf=[35 36 37 38];
rfc=[5 6 7];
rc=[12 13 14];
rcp=[19 20 21];
rp=[52 53 54 55];
rpo=[50 60 63];
rt=[40 42 44 46];

areas={laf;lf;lfc;lc;lcp;lp;lpo;lt;raf;rf;rfc;rc;rcp;rp;rpo;rt};

function [thSubjData]=thresholdData(subjData,density)
[nSubjs,d2,d3]=size(subjData);
thSubjData=zeros(nSubjs,d2,d3);
for i=1:nSubjs
    m=squeeze(subjData(i,:,:));
    thM=threshold_proportional(m,density);
    %thM(thM>0)=1;
    thSubjData(i,:,:)=thM;
end